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m 

O 

o 

(N 



Nikolay Prokof'ev and Boris Svistunov 
Department of Physics, University of Massachusetts, Amherst, MA 01003, USA and 
Russian Research Center "Kurchatov Institute", 123182 Moscow, Russia 

We report results of large-scale Monte Carlo simulations of superfluid-insulator transitions in 
commensurate 2D bosonic systems. In the case of off-diagonal disorder, we find that the transition 
is to a gapless incompressible insulator, and its dynamical critical exponent is z = 1.65 ±0.2. In the 
case of diagonal disorder, we prove the conjecture that rare statistical fluctuations are inseparable 
from critical fluctuations on the largest scales and ultimately result in the crossover to the generic 
universality class (apparently with z = 2). However, even at strong disorder, the universal behavior 
sets in only at very large space-time distances. This explains why previous studies of smaller clusters 
mimicked a direct superfluid-Mott-insulator transition. 
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Quantum phase transitions in disordered systems re- 
main a poorly understood phenomenon despite enormous 
interest in this field. The T — transition between the 
superfluid (SF) and insulating (I) phases is believed to 
determine properties of various condensed matter sys- 
tems: '^He in porous media and aerogels P, 0, ''He on 
various substrates P, 0, Q , thin superconducting films 
[slIR ItLIsL 1^ ITol lllll ■ Josephson-junction arrays |0, dis- 
ordered magnets etc. 

There are strong arguments that the basic Hamilto- 
nian which captures all the essential physics of the SF-I 
transition is the bosonic Hubbard model with disordered 
chemical potential [HITI IitIITsI I . In the limit of large 
occupation numbers, the bosonic Hubbard Hamiltonian 
is equivalent to the system of coupled Josephson junc- 
tions. Fermionic systems map to this Hamiltonian under 
the assumption that Cooper pairs are preformed at fi- 
nite temperature, and the transition is driven only by 
quantum fiuctuations of the phase of the complex order 
parameter. To deal with granular superconductors one 
may also introduce disorder to hopping amplitudes. 

It was suggested in Ref. 16] that one has to consider 
only two competing insulating phases — the incompress- 
ible (gaped) Mott-insulator phase (MI), and the com- 
pressible gapless Bose glass (BG) phase. However, more 
recently it was argued that apart from the BG phase 
characterized as a compressible insulator with variable- 
range-hopping conductivity at finite temperature [T^ . 
there may exist other phases such as a Bose metal with 
finite conductivity in the T — > limit ^19] and an in- 
compressible Mott glass with the conductivity pseudogap 
[20|. Theoretical calculations for the strongly coupled 
SF-I critical point are notoriously difficult and are not 
based on well controlled approximations since localiza- 
tion and interaction effects cannot be separated |^ . 
Thus even the qualitative understanding of the phase di- 
agram is still under debate. In particular, it was argued 
in El El El m that MI and SF phases are al- 
ways separated by the BG phase at any finite disorder. 
However experiments 0|, most Monte Carlo simulations 



|25l|2a|27l, and other theories [23, |2g|, |30| present evi- 
dence in favor of a direct transition between MI and SF 
phases (in the case of commensurate filling of the lattice 
and not so strong disorder). In d = 1 this contradiction 
was apparently resolved using arguments based on rare 
statistical fiuctuations ,2^] (Lifshitz-Griffiths-McCoy sin- 
gularities 0,11113^), renormalization- group e qua tions 
[SJillll, and quantum Monte Carlo simulations [Sql. 

In this Letter, we numerically address the problem of 
the SF-I transition in a disordered commensurate 2D sys- 
tem. Our large-scale simulations based on the classical 
Worm Algorithm [s^l demonstrate the absence of the di- 
rect SF-MI transition. We clearly see, however, that — 
even at strong disorder — the universal asymptotic long- 
range behavior sets in only at large space-time distances. 
This result, on one hand, explains previous observations 
of the direct SF-MI transitions in simulations of much 
smaller clusters, and, on the other hand, implies that 
the superfluid stiffness and compressibility should obey 
generic scaling laws only in a very close vicinity of the 
phase transition point which may be hard to study ex- 
perimentally. 

The Worm algorithm (WA) is a high-performance 
universal Monte Carlo scheme applicable to any model 
with the configuration space of continuous paths. The 
principle of WA is to work in an extended configuration 
space combining the physical closed-path sector and the 
broken-path sector with two path endpoints. All up- 
dates in the extended configuration space are through 
the motion of the endpoints (or even just one of them), 
so that the configuration evolution looks like a creep- 
ing worm. Though the WA updating strategy is based 
on local Metropolis moves (which is a key issue for its 
universality) , it has a remarkable efficiency equivalent to 
that of the Ijest cluster algorithms |^ . 

WA was introduced initially for the continuous-time 
path integral, or worldline, representation of lattice quan- 
tum systems [s^ . More recently it was implemented 
within the Stochastic Series Expansion (SSE) method 
[3^ . and generalized to classical lattice systems in the 
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closed-path representations jSgl- The optimal choice of 
WA depends on the problem being addressed; e.g. the 
quantum worldline scheme is best suited for ab initio 
simulations of bosons in optical lattices while the 
SSE scheme has certain storage and computational ad- 
vantages when simulating quantum models with a re- 
stricted Hilbert space, like spins and hard-core bosons 
[41I [4^. However, if one is interested in generic prop- 
erties of quantum phase transitions, then the optimal 
choice is a (c? -f l)-dimensional classical scheme, which 
is algorithmically superior from all points of view. This 
approach was advocated in Ref. 18], and most recently, 
using WA, in Ref. j4j|. 

The easiest way to derive the classical {d + 1)- 
dimensional equivalent of the d-dimensional quantum 
model is to start from the lattice path integral for par- 
ticle trajectories in imaginary time. The basic step is to 
"roughen" trajectories by replacing integrals over kink 
positions (particle jumps between lattice sites) in time by 
discrete sums. To this end one introduces a grid of imag- 
inary times and requires that kinks occur only at time 
slices forming the grid. Bonds in the spatial direction of 
thus obtained {d + l)-dimensional lattice carry an inte- 
ger charge — a spatial current equal to an algebraic sum of 
kinks associated with a given pair of sites and a time slice 
(the kink sign specifies its direction). The finite value of 
the time interval between slices, Ar, makes the roughen- 
ing procedure, generally speaking, not unique. Only in 
the physical limit of At the absolute values of spa- 
tial currents are either zero or unity. At strong roughen- 
ing the constraint that the absolute value of the spatial 
bond charge is < 1 is no longer necessary. The maximal 
value now depends on how we roughen the original tra- 
jectories with more than one kink between a given pair 
of sites during the time interval Ar — either as irrelevant 
rare events that are neglected, or by ascribing all such 
kinks to the same spatial bond of the {d+ l)-dimensional 
lattice. Note, there is no qualitative difference between 
the two cases, and the choice is just a matter of con- 
venience. Bonds in the temporal direction also carry an 
integer charge, equal to the occupation number of a given 
site between adjacent time slices. For the sake of symme- 
try, it is convenient to refer to the temporal-bond charges 
as temporal currents, so that each bond of the {d + 1)- 
lattice carries an integer current, J. The conservation 
of the particle number imposes an obvious constraint on 
bond currents — the divergence of the (d -I- l)-current at 
any {d + l)-site is zero. Graphically, if one represents 
the bond currents by oriented lines (in accordance with 
the sign of J) then all contributions to the partition- 
function will have a form of closed-path configurations 
of such lines, while configurations for the Green-function 
will have two endpoints. WA simulates both quantities 
simultaneously, by switchin g b etween partition-function 
and Green- function sectors [33, ■ 

Let us denote by { Jx,a} the bond current configuration 
where x = (r, r) are discrete space-time coordinates, and 
index a = fi, . . . , f^, f stands for unit vectors of axis di- 



rections, so that (x, a) defines a bond in the direction 
a, adjacent to the site x. As usual, the configuration 
weight W^[{Jx,a}] may be formally written as a Gibbs 
factor exp{— iJ/T} which for positive-definite W defines 
the classical bond Hamiltonian H/T = — In W. For mod- 
els with on-site particle-particle interactions the config- 
uration weight is simply given by the product of bond 
weights Wra{Jx,a), and, correspondingly. 
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The zero-divergence constraint can be written as 
J2a J^.a + J2a J^.-a = 0, whcrc, by definition, the di- 
rection —a is understood as opposite to a and J-x.,-a = 

In this paper, we are interested only in the universal 
critical behavior of the model JQ), which is supposed to 
be insensitive to quantitative details of the Hamiltonian. 
This freedom may be used to make spatial and tempo- 
ral directions symmetric with respect to each other. In 
the original quantum model, the temporal direction is 
not symmetric even with respect to its opposite because 
occupation numbers are always positive. Nevertheless, 
one may count them from some integer tiq, substitute 
>/x,f Jx,f — no, and formally consider Jx,f £ (^cOi oo)- 
Alternatively, one may introduce a symmetric constraint, 
say, |J| < 1. The latter case, reminiscent of the map- 
ping between bosonic and spin-1 systems, seems to be 
more natural and computationally economic. Histori- 
cally, however, a model with J G (—00, 00), motivated by 
its derivation from the Josephson-junction array Hamil- 
tonian fl8i] . was adopted. To retain the possibility of 
a direct comparison with previous numeric studies, we 
work with the same model 181: 
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In terms of the underlying bosonic system, K repre- 
sents the particle hopping amplitude in units of the on- 
site repulsion, the discrete field /ir = Mo + Mr involves 
the chemical potential mo and the white-noise diago- 
nal disorder fir G [—A, A]. In this study we are con- 
cerned with the commensurate filling of the lattice, i.e. 
n — {{Jx.f)) =integer, where ((. . .)) stands for the aver- 
age over all lattice points, statistical and disorder fluctu- 
ations, and thus set /io = 0. [An accurate study of the 
half-integer n case has been reported recently by Alet and 
S0rensen ^3]. We also consider model l(2Jl with the off- 
diagonal disorder introduced by letting the parameter K 
to be dependent on r and spatial direction, and confine 
ourselves to the case of a broken-bond disorder, where 
for some randomly chosen r and a' — fi,...fd we set 
K-rra' ~^ (equivalent to a rigid constraint J = on the 
corresponding bond). The phase diagram of the homo- 
geneous system is shown in Fig.^ the SF-I transition at 
the commensurate filling is located at K^^ — 0.33305(5) 
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FIG. 1: The phase diagram of model Q in the absence of 
disorder. Error bars are of order 10~^ and smaUer than point 
size. 
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FIG. 2: Correlation functions G(r, 0) and G(0, r) for K = 
0.3810, 0.3813, 0.3815 and system size L x L x = 160 x 
160 X 500. Error bars are comparable to the line widths. 



[43 ■ The value of the MI gap is defined as half the differ- 
ence between the critical values of the chemical potential 
in Fig. in i.e. 2i?gap(if) = ^i^P^ - /ii'*"""'. 

We start with the off-diagonal disorder case, and con- 
sider a system with a quarter of all bonds being broken. 
Typically, we include about 10'^ disorder realizations into 
the statistics for system sizes L < 40, and 4 x 10'*/i 
for larger L. The critical point, K^-, and the dynami- 
cal exponent, z, may be obtained from the study of the 
Green function, G(r,T), naturally evaluated within the 
WA approach [33|. At the critical point one should see 
a power-law decay: G{r,0) r^i^+v) as r ^ 00, and 
G(0,t) T--(i+»;/2^) as T ^ 00. This way we find (see 
Fig. 121 as weh as Figs. 01 andO 



FIG. 3: Suporfluid stiffness of the broken-bond model as a 
function of K at different system sizes; 40 x 40 x 40 - open 
circles, 80 x 80 x 80 - filled circles, 160 x 160 x 160 - open 
squares, 160 x 160 x 500 - filled squares. 
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FIG. 4: Compressibility of the broken-bond model as a func- 
tion of K at different system sizes; 40 x 40 x 40 - open circles, 
80 X 80 X 80 - filled circles, 160 x 160 x 160 - open squares, 
160 X 160 X 500 - filled squares. 



z + 7] = 1.37(8), 1 + = 0.82(6), i.e. 

z = 1.65(20) , 
77 = -0.3(1) . 



(4) 
(5) 



K,, = 0.3813(2) , 



(3) 



It is clear in Fig. [3 that the asymptotic behavior of 
the correlation function sets in only at sufficiently large 
space-time distances > 10 lattice periods. Moreover, the 
short-range behavior of G mimics the critical point of 
the SF-MI transition in the regular system, where z = 1. 
This peculiar behavior implies that the curves for the su- 
perfluid stiffness, p^, and compressibility, k, will acquire 
their universal forms only in a very narrow region around 
the critical point. Away from this region, the form of 
Ps(if — and k,(K — K^) curves should be essentially 
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FIG. 5: Superfluid stiffness for the diagonal disorder case as 
a function of K at different system sizes; 10 x 10 x 20 - filled 
circles, 20 x 20 x 49 - open circles, 40 x 40 x 121 - filled squares, 
80 X 80 X 298 - open squares, 160 x 160 x 733 - triangle down 



different, as suggested by the extended transient evolu- 
tion of z from « 1 to its true critical value. In Figs. |21 and 
0]we indeed observe such a behavior. The anomalously 
narrow critical region makes it virtually impossible — even 
with our large cluster sizes — to reliably determine the 
correlation radius critical exponent v. Along with the dy- 
namical exponent z it is supposed to determine the crit- 
ical behavior of the compressibility, k oc {K — KcY^^^^^f 
and superfluid stiffness, ps oc {K - KcY^ |]Jj. The data 
in Figs. El and 0] at best guarantee only the inequalities 
v{2 — z) < 1 and vz > 1, but do not allow us to test the 
Harris criterion 0| v > 2/d = 1. 

The finite-size scaling of the data for compressibility 
demonstrates no sign of saturation below Kc and thus 
strongly suggests that in the insulating state the com- 
pressibility vanishes. Though the insulating state is in- 
compressible, it is easy to prove that it is gapless and 
thus is qualitatively different from the conventional Mott 
insulator and Bose Glass states. Indeed, in an infinite 
system it is always possible to find an arbitrarily large 
cluster that is nearly uniform (in the sense that fluctua- 
tions of K away from its cluster average value are arbi- 
trarily small/rare). Taking into account that Kc in the 
disordered system is larger than the ideal-system critical 

value k'"c'\ we conclude that such clusters are nothing 
else but finite-size superfluid lakes. Hence, the gap asso- 
ciated with adding one more particle to the cluster scales 
as l/l'^, where / is the cluster size. The absence of an 
upper bound on I immediately implies the absence of the 
global gap in the system spectrum and finite optical con- 
ductivity. 

The diagonal-disorder case is different. Previously re- 
ported data j23| for small clusters L x L < 12 x 12 and 
disorder strength A = 0.2 were interpreted as a direct 
SF-Ml transition with z = 1. We extended the study 



FIG. 6: Compressibility for the diagonal disorder case as a 
function of K at different system sizes; 10 x 10 x 20 - filled 
circles, 20 x 20 x 49 - open circles, 40 x 40 x 121 - filled squares, 
80 X 80 X 298 - open squares, 160 x 160 x 733 - triangle down, 
160 X 160 X 160 - triangle up. The data for L > 80 collapse 
on each other within the error bars. 



of the A — 0.2 case to system sizes L x L < 160 x 160 
and did not find any deviations from the direct transi- 
tion picture with vanishingly small compressibility be- 
low ^^(A = 0.2) = 0.325(1). However, the value of 
the Ml gap in the ideal system is almost three times 
smaller than A aX K = Kc, see Fig. According to 
the argument/theorem of Refs. [l^ l24j. the state with 
A > Egap is a compressible (gapless) insulator, or EG, 
because in the infinite system one can always find ar- 
bitrary large regions with the chemical potential being 
nearly homogeneously shifted downwards or upwards by 
A (and thus they are doped with particles or holes). 
Since the distance between such regions is exponentially 
large for A ^ 0, their effect is simply undetectably small 
for A = 0.2. 

Even if the state right below Kc is a compressible insu- 
lator, the question remains whether Griffiths-McCoy sin- 
gularities are inseparable from critical fluctuations and 
ultimately result in the crossover to the generic SF-BG 
transition, or they merely provide a regular background 
contribution to k on which a singular contribution Kging 
is superimposed. The latter scenario implies a cusp on 
the compressibility curve, and criticality different from 
SF-BG. To answer this question we performed simula- 
tions for disorder strength A = 0.4. As before, the ideal 
Ml gap at the transition point Kc = 0.2910(5) is about 
two times smaller than A, and k has to be finite at Kc- 

In Figs.|31and|Slwe show the data for the compressibil- 
ity and superfluid stiffness which away from the critical 
point mimic the ideal-system behavior (with the corre- 
lation length exponent v ~ 0.7 and z « 1), but close to 
Kc show a spectacular crossover to another universality 
class. Strong finite-size corrections to k for system sizes 
L < 20 saturate for L > 20, and the thermodynamic 
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FIG. 7: A map of doped places slightly below Kc showing 
the picture of rare, well isolated regions. Point sizes are pro- 
portional to <j)M(jc). 

curve clearly demonstrates finite^ and non-singular de- 
pendence k{K — Kc). At the same time, we observe a 
crossover in the ps{k — Kc) dependence, and see that ps 
approaches zero with zero derivative, i.e. vz > 1. From 
the decay of the Green function at the critical point we 
obtain 

z = 2.0(2) , (6) 
V = 0.11(2) . (7) 

Unfortunately, the large-scale crossover did not allow us 
to determine the critical exponent from this set of data. 
Recent data for half-integer n are best fit with v = 1.15, 
but they also suffer from large finite-size corrections ^3 ■ 
Apparently, the best strategy in the future is to search 
for a classical model with the smallest crossover scale. 

Within the WA approach one may directly visualize 
places where particles are added/removed at low tem- 
perature in a sample with a given disorder realization. 
The standard procedure of subtracting density maps ob- 
tained in canonical simulations n{i,N) — n{i, N — I) is 
time consuming, and for large systems requires extremely 
high-precision data for n(i,N). The new technique is 
based on the statistics of the Green function calculated 
at the chemical potential between the steps on the N{p) 
curve, i.e. at ^ = Eg{N) - Eg{N - 1), where Eg{N) is 
the ground state energy of the A^-particle system. In the 
T ^ limit we write G{N;r,r',T = 13/2) = F{N;r,r') 
and notice that only ground states with — 1 and N 
particles contribute to the answer. Thus F{N;r,r') w 
(/>jv(i")0Ar(r') where 0Ar(r) is given by the positive-definite 
groundstate-groundstate matrix element 

^N{r)^{^G{N)\bl\'^G{N ~1)) , (8) 

and may be viewed as the added-particle "wavefunction" . 
Its localization length diverges at the SF-I transition. 



The two important parameters which characterize the 
structure of the normalized wavefunction are the local- 
ization radius 

n'^{[r-{r)f)=J2 [r-(r)]^0^(r), (9) 

r 

and the state "area" , or the number of sites over which 
(f) delocalizes, 

The dependence of ^ on 7^ gives the fractal dimension of 
the state. It is also important to monitor correlations in 
the overlaps between different states, (/)Ar(r)0Ar/ (r). 
For example, the percolation type scenario of the SF- 
I transition assumes large fractal superfluid "lakes" in 
the insulating phase; if so, then with the probability of 
order unity there must exist an almost complete overlap 
between a pair of states (</'jv„ , ^jVi,) with Na, from a 
narrow interval of width SN <^ N around N. 

In Fig.[7|wc show a typical map of J^r' = 1; r, r') 
for the insulating state a.t K — 0.288 and system size 
L X L X Lr — 160 X 160 x 1000. One can clearly see 
isolated regions doped with particles. 

In summary, we have performed large-scale simula- 
tions of the superfluid-insulator transition in the (2-1-1)- 
dimensional classical analog of the commensurate disor- 
dered 2D bosonic system. For diagonal disorder, our re- 
sults suggest that commensurability is not relevant in the 
long-range limit, and the universality class of the tran- 
sition (superfluid-Bose-glass) is the same for all filling 
factors. In particular, we unambiguously resolved the fi- 
nite compressibility at the critical point originating from 
rare statistical fluctuations and demonstrated that they 
are responsible for the crossover in critical behavior. We 
believe that our results rule out the earlier-reported di- 
rect superfiuid-Mott-insulator transition in this model. 

In the off-diagonal disorder case, the compressibility 
vanishes at the critical point. The incompressible in- 
sulating phase, however, is gapless, and its universality 
class is characterized by the dynamical critical exponent 
z = 1.65 ±0.2. 

A general observation is that even for large diagonal 
and off-diagonal disorder, the universal asymptotic long- 
range behavior sets in only at large space-time distances 
{r-^ 20 lattice periods) . This circumstance explains previ- 
ous observations of the direct superfluid-Mott-insulator 
transition in small-size clusters and implies that the 
the superfluid stiffness and compressibility should obey 
generic scaling laws only in a very close vicinity of the 
phase transition point. 
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